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Abstract 

The purpose of this article is to provide a starter kit for multicanonical simulations in statistical physics. Fortran 
code for the g-state Potts model in d = 2, 3, . . . dimensions can be downloaded from the Web and this paper 
describes simulation results, which are in all details reproducible by running prepared programs. To allow for 
comparison with exact results, the internal energy, the specific heat, the free energy and the entropy are calculated 
for the d — 2 Ising (g = 2) and the q — 10 Potts model. Analysis programs, relying on an all-log jackknife technique, 
which is suitable for handling sums of very large numbers, are introduced to calculate our final estimators. 
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1. Introduction 

The conventional Metropolis [1] method simu- 
lates the Gibbs canonical ensemble at a fixed tem- 
perature T and allows for easy calculations of the 
(internal) energy and functions thereof. However, 
some of the most important quantities of statisti- 
cal physics, free energy and entropy, can only be 
obtained by tedious integrations. One way to over- 
come this problem is by multicanonical (MUCA) 
simulations, which calculate canonical expectation 
values over a temperature range in a single simu- 
lation by using the weight factor [2] 
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where n(E) is the number of states with energy E 
and S(E) the microcanonical entropy. In an exten- 
sion of the microcanonical terminology on may call 
b(E) microcanonical inverse temperature (b(E) = 



1/T(E) in natural units with Boltzmann constant 
ks = 1) and a(E) microcanonical, dimensionless 
free energy, see appendix A. 

MUCA simulations became popular with the in- 
terface tension calculation [2] of the 2d 10-state 
Potts model, when the method emerged as the 
winner of largely disagreeing estimates, which af- 
ter their publication became resolved by exact val- 
ues [3] . Similar simulation concepts can actually be 
traced back to the work by Torrie and Valleau [4] 
in the 1970s. In recent years the MUCA method 
has found many applications, besides for first or- 
der phase transitions mainly for complex systems 
including spin glasses and peptides, see [5] for a 
brief review and a summary of related methods. 

The scope of this article is limited to the Ising 
model and its generalization in form of q-state 
Potts models, for a review see [6]. Fortran rou- 
tines which work in arbitrary integer dimensions 
d = 2, 3, . . . are provided, but we confine our 
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demonstrations to d = 2, to allow for comparison 
with rigorous analytical calculations. The Ising 
model simulation is seen to match the exact finite 
lattice results of Ferdinand and Fisher [7], while 
for the q = 10 Potts model one finds agreement 
with the rigorously known transition temperature 
and latent heat of Baxter [8] . Details of the model 
are summarized in the first part of section 2. In 
the second part of this section the downloading of 
the Fortran code and its use are explained. 

In section 3 MUCA simulations are treated. The 
temperature dependence of the standard thermo- 
dynamic quantities - energy specific heat, free en- 
ergy and entropy - is calculated for the 2d Ising 
model as well as for the 2d 10-state Potts model and 
the canonically re- weighted histograms are shown. 
Special attention is given to the analysis procedure, 
which has to be able to handle sums of very large 
numbers. This is done by using only the logarithms 
until finally the quotient of two such numbers is 
obtained. Jackknife [9] binning is used to minimize 
bias problems which occur in the re-weighting of 
the simulation data to canonical ensembles. 

Using the provided Fortran code and following 
the instructions allows for a step by step reproduc- 
tion of the figures and all other numerical results 
presented in this article. This could be a desirable 
standard for more involved simulations too. Some 
conclusions are given in the final section 4. 



2. Getting Started 

2.1. The Potts model 

We introduce the Potts models on <i-dimcnsional 
hypercubic lattices with periodic boundary condi- 
tions. For this paper we stay close to the notation 
used in the accompanying computer programs and 
define the energy E via the action variable 
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where 5(qi, qj) is the Kronecker delta. The sum < 
ij > is over the nearest neighbor lattice sites and 
q^ is the Potts state of configuration k at site i. 



(k) 

For the g-state Potts model q\ ' takes on the val- 
ues 1, . . . , q. As the variable iact takes on integer 
values, it allows for convenient histograming of its 
values during the updating process. Occasionally, 
we use the related mean values 



actm = iact / (dN) and e s = E / N . 



(3) 



Each configuration (microstate of the system) k 
defines a particular arrangements of all states at 
the sites and, vice versa, each arrangement of the 
states at the sites determines uniquely a configu- 
ration: 
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The expectation value of an observable O is defined 
by 
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where the sum is over all microstates and the par- 
tition function Z — Z(0) normalizes the expecta- 
tion value of the unit operator to (1) = 1. As there 
are q possible Potts states at each site, the total 
number of microstates is 



Z(/3 = 0) = q 



N 



(6) 



Including = in a MUCA simulation allows for 
the normalization of the partition function neces- 
sary to calculate the canonical free energy and the 
entropy as a function of the temperature. 

Our definition(5) of agrees with the one com- 
monly used for the Ising model [10], but disagrees 
by a factor of two with the one used for the Potts 
model in [8,3]: 



= /3 Ising = Potts /2. 
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For the 2d Potts models a number of exact results 
are known in the infinite volume limit. The critical 
temperature [8] is 

I = & = = \ ln(l + ^), g = 2,3,... .(8) 

The phase transition is second order for q < 4 and 
first order for q > 5. At C the average energy per 
Potts state is [8] 



-e c s = 2 + 2/Jq 



(9) 
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Fig. 1. Fortran code directory structure. 

where, by reasons of consistency with the Ising 
model notation, also our definition (3) of e s differs 
by factor of two from the one used in most Potts 
model literature. For the first order transitions at 
q > 5 equation (9) gives the average of the limit- 
ing energies from the ordered and the disordered 
phase. The exact infinite volume latent heat Ae s 
and the entropy jumps As were also calculated by 
Baxter [8] , whereas the interfacial tensions f s were 
derived more recently [3] . 

2.2. The Fortran Code 

Figure 1 shows the directory tree in which the 
Fortran routines are stored. MUCA is the parent di- 
rectory and on the first level we have the directo- 
ries Exercises, ForLib, ForProg and Work. The 
master code is provided in the directories ForLib 
and ForProg. ForLib contains the source code of 
a library of functions and subroutines. The library 
is closed in the sense that no reference to non- 
standard functions or subroutines outside the li- 
brary is ever made. The master versions of the main 
programs and certain routines (which need input 
from the parameter files discussed below) are con- 
tained in the subdirectory ForProg. The demon- 
strations of this article are contained in the subdi- 
rectories of Exercises. 

To download the code, start with the URL 
www.hep.fsu.edu/~berg and click the Research 
link, then the link Mult i canonical Simulations. 
On this page follow the link Fortran Code and 
get either the file muca.tar (399KB) or the file 
muca.tgz (40KB). On most Unix platforms you 
obtain the directory structure of figure 1 from 
muca.tar by typing 



alternatively from muca.tgz by typing either tar 
-zxvf muca.tgz or gunzip muca.tgz followed by 
(10). 

You obtain the results of this paper by compil- 
ing and running the code prepared in subdirecto- 
ries of Exercises, e.g. f77 -0 program, f followed 
by ./a. out. Due to the include and parameter 
file structure used, the programs and associated 
routines of ForProg compile only in subdirectories 
which are two levels down from the MUCA parent 
directory. This organization should be kept, unless 
you have strong reasons to change the dependen- 
cies. The present structure allows to create Work di- 
rectories for various projects, with the actual runs 
done in the Work subdirectories. Note that under 
MS Windows with the (no longer marketed) MS 
Fortran compiler a modification of the include 
structure of the code turned out to be necessary. If 
such problems are encountered, one solution is to 
copy all needed files to the subdirectory in question 
and to modify all include statements accordingly. 

Each Exercises subdirectory contains a 



readme.txt 



(11) 



file with instructions, which should be followed. 
The subdirectories el. . . prepare some code to 
check for the correctness of the conventional 
Metropolis code and the subdirectories e2. . . 
prepare examples of multicanonical simulations, 
which are discussed in the next section. The sim- 
ulation parameters are set in the files 

lat.par, lat.dat, potts. par, mc.par and muca.par ,(12) 

or a subset thereof, which are also kept in each of 
the subdirectories of Exercises. The dimension 
nd of the system and the maximum lattice length 
ml are defined in lat .par. In lat .dat the lattice 
lengths for all directions are assigned to the ar- 
ray nla, which is of dimension nd. This allows for 
asymmetric lattices. The number of Potts states, 
nq, is defined in potts. par. The parameters of 
the conventional Metropolis simulation are defined 
in mc .par: These are the (3 value beta, which de- 
fines the initial weights in case of a MUCA sim- 
ulation, the number of equilibrium sweeps nequi, 
the number of measurement repetitions nrpt and 
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the number of measurement sweeps nmeas. Addi- 
tional parameters of the MUCA recursion are de- 
fined in muca.par: The maximum number of re- 
cursions nrecjnax, the number of sweeps between 
recursions nmucasw and the maximum number of 
tunnelings (16) maxtun, which terminates the re- 
cursion unless nrecjnax is reached first. 

Whenever data for a graphical presentation are 
generated by our code, it is in a form suitable for 
gnu-plot, which is freely available. Gnuplot driver 
files fin. pit are provides in the solution directo- 
ries, such that one obtains the the plot by typing 
gnuplot fin. pit on Unix and Linux platforms 
(under MS Windows follow the gnuplot menu) . 



3. Multicanonical Simulations 

A conventional, canonical simulation calculates 
expectation values at a fixed temperature T and 
can, by re- weighting techniques, only be extrapo- 
lated to a vicinity of this temperature [11]. In con- 
trast, a single MUCA simulation allows to obtain 
equilibrium properties of the Gibbs ensemble over 
a range of temperatures, which would require many 
canonical simulations. This coined the name multi- 
canonical. The MUCA method requires two steps: 

(i) Obtain a working estimate w mu (k) of the 
weights W\/ n {k). Working estimate means 
that the approximation to (1) has to be good 
enough to ensure movement in the desired 
energy range, but deviations of w mu (E) from 
(1) by a factor of, say, ten are tolerable. 

(ii) Perform a Markov chain MC simulation with 
the final, fixed weights w mu (k). Canonical 
expectation values are found by re- weighting 
to the Gibbs ensemble. 

To obtain working estimates w mu (k) of the 
weight factors (1), a slightly modified version of 
the recursion of Rcf. [12] is used. As the analytical 
derivation of the modified recursion has so far only 
been published in conference proceedings, it is for 
the sake of completeness included in appendix A 
of this paper. The Fortran implementation is given 
by the subroutine 

p_mu_rec.f (13) 



of ForProg. One subtlety is that two histogram 
arrays hup and hdn are introduced to keep sepa- 
rately track of the use of upper and lower entries 
of nearest neighbor pairs. Detailed explanations of 
the code will be part of a book [13]. 

The question whether more efficient recursions 
exists is far from being settled. For instance, F. 
Wang and Landau [14] made recently an interest- 
ing proposal. Exploratory comparisons with the 
recursion used in the present paper reveal similar 
efficiencies [15]. 

In between the recursion steps the Metropolis 
updating routine 

potts_met_f (14) 

is called, which implements the standard Metropo- 
lis algorithm [1] for general weights. The random 
number generator of Marsaglia et al.[16] is inte- 
grated to ensure identical results on distinct plat- 
forms. 

3.0.1. Example runs 

First, we illustrate the MUCA recursion for the 
20 2 Ising model. We run the recursion in the range 

namin = 400 < iact < 800 = namax . (15) 

These values of namin and namax are chosen 
to cover the entire range of temperatures, from 
the completely disordered (fi — 0) region to the 
groundstate (/3 — > oo). In many applications the 
actual range of physical interest is smaller, namin 
and namax should correspondingly be adjusted, 
because the recursion time increases quickly with 
the range. The recursion is completed after maxtun 
tunneling events have been performed. A tunnel- 
ing event is defined as an updating process which 
finds its way from 

iact = namin to iact = namax and back .(16) 

This notation comes from the applications of the 
method to first order phase transitions [2], for 
which namin and namax are separated by a free en- 
ergy barrier in the canonical ensemble. Although 
the MUCA method removes this barrier, the ter- 
minus tunneling was kept. The requirement that 
the process tunnels also back is included in the 
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definition, because a one way tunneling is not in- 
dicative for the convergence of the recursion. Most 
important, the process has still to tunnel when 
the weights are frozen for the second stage of the 
simulations Note that things work differently for 
the Wang-Landau recursion [14]. It has no prob- 
lems to tunnel in its initial stage, but its estimates 
of the spectral density are still bad, such that the 
tunneling process gets stuck as soon as the weights 
are fixed. 

For most applications ten tunnelings during our 
recursion part lead to acceptable weights. If the 
requested number of tunnelings is not reached af- 
ter a certain maximum number of recursion steps, 
the problem will disappear in most cases by rerun- 
ning (eventually several times) with different ran- 
dom numbers. Otherwise, the number of sweeps be- 
tween recursions should be enlarged, because our 
recursion is strictly only valid when the system is 
in equilibrium. One may even consider to discard 
some sweeps after each recursion step to reach equi- 
librium, but empirical evidence indicates that the 
improvement (if any) does not warrant the addi- 
tional CPU time. The disturbance of the equilib- 
rium is weak when the weight function approaches 
its fixed point. In the default setting of our pro- 
grams we take the number of sweeps between recur- 
sion steps inversely proportional to the acceptance 
rate, because equal number of accepted moves is a 
better relaxation criterion than an equal number 
of sweeps. 

In the subdirectory e2_01 of Exercises a 20 x 
20 lattice Ising model simulation is prepared for 
which we requested ten tunneling events. We we 
find them after 787 recursions and 64,138 sweeps, 
corresponding to an average acceptance rate of 
20*787/64138=0.245 (the acceptance rate can be 
calculated this way, because the number of ac- 
cepted sweeps triggers the recursion). Almost half 
of the sweeps are spent to achieve the first tunnel- 
ing event. Subsequently, an MUCA production run 
of 10,000 equilibrium and 32 x 10, 000 sweeps with 
measurements is carried out. On a GHz Linux PC 
the entire runtime (recursion plus production) is 
about thirty seconds. In the subdirectory e2_02 a 
similar simulation is prepared for the 2d 10-state 
Potts model on a 20 2 lattice. 



3.1. Re-weighting to the canonical ensemble 

Let us assume that we have performed a MUCA 
simulation which covers the action histogram 
needed for a temperature range 

flnin < (i < /3max • (17) 

In practice this means that the parameters namax 
and namin in muca . par have to be chosen such that 

namin <C act(/3 m ; n ) and act(/3 max ) <C namax 

holds, where act(/3) is the canonical expectation 
value of the action variable (2) . The <C conditions 
may be relaxed to equal signs, if b(iact) = /3 m ; n 
is used for all action values iact < act(/3 m i n ) 
and 6(iact) = /3 max for all action values iact > 
act(/3 max ). 

Given the MUCA time series, where i = 1, . . . , n 
labels the generated configurations, the defini- 
tion (5) of the canonical expectation values leads 
to the MUCA estimator 

0= (18) 

£? =1 0« exp [-0E® + 6(£?W) £?(0 - a (£«)] 

Er=i ex P [-P E(i) + b(E(*)) EM - a(E(*))] ■ 

This formula replaces the MUCA weighting of 
the simulation by the Boltzmann factor of equa- 
tion (5). The denominator differs from Z by a con- 
stant factor, which drops out because the numer- 
ator differs by the same constant factor from the 
numerator of (5). If only functions of the energy 
(in our computer programs the action variable) 
are calculated, it is sufficient to keep histograms 
instead of the entire time series. For an operator 
OW = /(£«) equation (18) simplifies then to 

7= (19) 

Ee /Cg) Ug) exp [-0 E + b(E) E - a(E)} 
Y jE h mu {E)e^[-0E + b{E)E-a{E)] 

where h mu (E) is the histogram sampled during 
the MUCA production run and the sums are over 
all energy values for which h mu (E) has entries. 
When calculating error bars for estimates from 
equations (18) or (19), we employ jackknife [9] es- 
timators to reduce bias problems. 
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A computer implementation of equations (18) 
and (19) requires care. The differences between the 
largest and the smallest numbers encountered in 
the exponents can be really large. To give one ex- 
ample, for the Ising model on a 100 x 100 lattice 
and (3 — 0.5 the groundstate configuration con- 
tributes —/3E = 10 4 , whereas for a disordered con- 
figuration E = is possible. Clearly, overflow dis- 
asters will result, if we ask Fortran to calculate 
numbers like exp(10 4 ). When the large terms in 
the numerator and denominator take on similar or- 
ders of magnitude, one can avoid them by subtract- 
ing a sufficiently large number in all exponents of 
the numerator as well as the denominator, result- 
ing in a common factor which divides out. Instead 
of overflows one encounters harmless underflows 
of the type exp(— 10 4 ). We implement the idea in 
a more general fashion, which remains valid when 
the magnitudes of the numerator and the denom- 
inator disagree. We avoid altogether to calculate 
large numbers and deal only with the logarithms 
of sums and partial sums. 

We first consider sums of positive numbers and 
discuss the straightforward generalization to arbi- 
trary signs afterwards. For C = A + B with A > 
and B > we calculate In C = ln(A + B) from the 
values In A and In B, without ever storing either A 
or B or C. The basic observation is that 



lnC = ln 



max(A, B) 1 + 



min(A, B) 



(20) 



Energy per spin e s 
Multicanonical data 



-0.5 



-1.5 



max(A, B) 

= max (In A, In B) + 
ln{ 1 + cxp [min(ln A, In B) — max(ln A, In B)\ } 

holds. By construction the argument of the expo- 
nential function is negative, such that an underflow 
occurs when the difference between min(ln A, In B) 
and max(ln^4, InB) becomes too big, whereas it 
becomes calculable when this difference is small 
enough. 

To handle alternating signs one needs in addition 
to equation (20) an equation for In |C| = In \A — B\ 
where A > and B > still holds. Assuming 
\nA ^ InS, equation (20) converts for ln|C| = 
In \A - B\ into 

In |C| = max (In A, In 5) (21) 
+ ln{ 1 - exp [min(ln A, In B) - max(ln A, In B)] } 



0.1 0.2 0.3 0.4 0.5 0.6 



Fig. 2. Energy per spin e s versus f3 for the 2d Ising model 
on an 20 X 20 lattice, Multicanonical data are compared 
with the exact result of Ferdinand and Fisher (full line), 
see the subdirectory e2_02. 

and, because the logarithm is a strictly monotone 
function, the sign of C = A — B is positive for 
In A > In B and negative for In A < In B. 

The computer implementation of equations 
(20) and (21) is provided by the Fortran function 
addln . f and the Fortran subroutine addln2 . f of 
ForLib, respectively. The subroutines potts_zln . f 
and potts_zln0.f of ForLib rely on this to per- 
form the jackknife re- weighting analysis for various 
physical observables. 

3.2. Energy and specific heat calculations 

We are now ready to analyze the MUCA data 
for the energy per spin of the 2d Ising model on a 
20 x 20 lattice, which we compare in figure 2 with 
the exact results of Ferdinand and Fisher [7] . The 
code is prepared in the subdirectory e2_03. 

The same numerical technique allows us to to 
calculate the specific heat, which is defined by 



C = 



dE 
If 



= (? ((E^)-(Ef) . 



(22) 



Figure 3 compares the thus obtained MUCA data 
with the exact results of Ferdinand and Fischer [7] . 
The code is also prepared in subdirectory e2_03. 

Figure 4 shows the energy histogram of the 
MUCA simulation together with its canonically 
re- weighted descendants at/3 = 0,/3 = 0.2 and 
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Fig. 3. Specific heat versus j3 for the 2d Ising model on 
a 20 X 20 lattice. Multicanonical data are compared with 
the exact result of Ferdinand and Fisher (full line), see the 
subdirectory e2_03. 




Fig. 4. Energy histogram from a multicanonical simulation 
of the 2d Ising model on a 20 X 20 lattice together with 
the canonically re-weighted histograms at /3 = 0, /3 = 0.2 
and (i = 0.4, see the subdirectory e2_03. 

f3 = 0.4. The Fortran code is prepared in the sub- 
directory e2_04. The normalization of the MUCA 
histogram is adjusted such that it fits reasonably 
well into one figure with the three re-weighted his- 
tograms. In figure 4 it is remarkable that the error 
bars of the canonically re- weighted histograms are 
not just the scale factors times the error bars of the 
MUCA histogram, but in fact much smaller. This 
can be traced to be an effect of the normalization. 
The sum of each canonical jackknifc histogram is 
normalized to the same number and this reduces 
the spread. 



Fig. 5. Multicanonical mean action variable (3) data for 
the 2d 10-state Potts model on a 20 X 20 lattice, see the 
subdirectory e2_02. 
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Fig. 6. Action histogram from a multicanonical simulation 
of the 2d 10-state Potts model on a 20 X 20 lattice together 
with the canonically re-weighted histograms at f3 = 0.71, 
see the subdirectory e2_06. 

Relying on the 2d 10-state Potts model data 
of the run prepared in subdirectory e2_02, we re- 
produce in subdirectory e2_05 the action variable 
actm results plotted in figure 5. Around (i = 0.71 
we observe a sharp increase of actm from 0.433 at 
(3 = 0.70 to 0.864 at f3 = 0.72, which signals the 
first order phase transition of the model. 

In figure 6 we plot the canonically re-weighted 
histogram at f3 = 0.71 together with the MUCA 
histogram using suitable normalizations, as pre- 
pared in subdirectory e2_06. The ordinate of fig- 
ure 6 is on a logarithmic scale and the canon- 
ically re-weighted histogram exhibits the double 
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Fig. 7. Free energies from multicanonical simulations of the 
2d Ising and 2d 10-state Potts models on a 20 X 20 lattice, 
see subdirectories e2_03 and e2_05. The lines are the exact 
results of Ferdinand and Fischer for the Ising model and 
the asymptotic equation (26) for the 10-state Potts model. 

peak structure which is characteristic for first order 
phase transitions. The MUCA method allows then 
to estimate the interface tension of the transition 
by calculating the minimum to maximum ratio on 
larger lattices, see [2,5]. 

3.3. Free energy and entropy calculations 

At f3 = the Potts partition function Z is 
given by equation (6). MUCA simulations allow 
for proper normalization of the partition func- 
tion by including f3 = in the temperature range 
(17). The normalized partition function yields im- 
portant quantities of the canonical ensemble, the 
Helmholtz free energy 



F = -/T 1 \n{Z) 
and the entropy 
S = P(F-E) 



(23) 



(24) 



where E is the internal energy (2). 

Figure 7 shows the free energy density per site 



/ = F/N 



(25) 



for the 2d Ising model as well as for the 2d 10-state 
Potts model. The Ising model analysis is prepared 
in the subdirectory e2_03 and the q = 10 analysis 
in e2_05. As in previous figures, the Ising model 



Fig. 8. Entropies from multicanonical simulations of the 
2d Ising and 2d 10-state Potts models on a 20 X 20 lattice, 
see the subdirectories e2_03 and e2_05. The full line is the 
exact result of Ferdinand and Fischer for the Ising model. 

data are presented together with the exact results, 
whereas we compare the 10-state Potts model data 
with the [3 — > oo asymptotic behavior. For large (3 
the partition function of our g-statc Potts models 
approaches q exp(—2f]dN+2(3dN/q) and, there- 
fore, 



/as 



2d 



2d- (3 



-i HQ) 



(26) 



q N 
Finally, in figure 8 we plot for the entropy density 



s = S/N 



(27) 



the 2d Ising model the MUCA results together with 
the exact curve. Further, entropy data for the 2d 
10-state Potts model are included in this figure. In 
that case we use s/3 instead of s, such that both 
graphs fit into the same figure. The analysis code 
for the entropy is contained in the same subdirec- 
tories as used for the free energy. 

In all the figures of this section excellent agree- 
ment between the numerical and the analytical re- 
sults is found. 



4. Conclusions 

There are many ways to extend the multicanoni- 
cal simulations of this paper. The interested reader 
is simply referred to the literature [5]. The pur- 
pose of this article is to serve a start-up kit for the 
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computer implementation of some of the relevant 
steps. There appears to be need for this, because 
to get the first program up and running appears 
to be a major stumbling block in the way of using 
the method. 

In the opinion of the author, multicanonical sim- 
ulations have the potential to replace canonical 
simulations as the method of first choice for stud- 
ies of small to medium-sized systems. As seen here, 
once the recursion necessary for the first part of a 
MUCA simulation is programmed, the entire ther- 
modynamics of the system follows from the second 
part of the simulation. However, the slowing down 
for larger system sizes is rather severe. 

Quite a number of similar methods exists, see [5] 
for a summary. A sound comparison would require 
that the goals of the simulations and their bench- 
marks arc defined first. So far the community has 
not set such standards. 
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Appendix A. Weight Recursion 

We first discuss the weights (1). By definition, 
the microcanonical temperature is 



b(E) 



1 



T(E) 



dS(E) 
dE 



(A.l) 



and we define the dimensionless, microcanonical 
free energy by 



a(E) = 



F{E) 



E 



T(E) T(E) 
= b(E)E - S(E) . 



-S(E) 



(A.2) 



It is determined by relation (A.l) up to an (irrele- 
vant) additive constant. We consider the case of a 
discrete minimal energy e and choose 



b(E) = [S(E + e) - S(E)} /e 
as the definition of b(E) . The identity 
S(E) = b(E) E — a(E) 



(A.3) 



implies 

S(E) - S(E - e) = 
b(E)E - b{E - e)(E - e) - a(E) + a(E - e) . 
Inserting eb{E - e) = S{E) - S(E - e) yields 

a(E - e) = a(E) + [b{E - e) - b(E)} E (A.4) 

and a(E) is fixed by defining a(E max ) = 0. Once 
b(E) is given, a(E) follows. 

A convenient starting condition for the initial 
(n = 0) simulation is 



b°(E) = and a°(E) = 0, 



(A.5) 



because the system moves freely in the disordered 
phase. Other b°(E) choices are of course possible. 
Our Fortran implementation allows for b° (E) — (3 
with [3 defined in the parameter file mc .par. 

The energy histogram of the n th simulation is 
given by H n (E). To avoid H n (E) = 0we replace 
for the moment 



H n {E) -» H n (E) = max [ho, H n (E)} , 



(A.6) 



where ho is a number < ho < 1. Our final equa- 
tions allow for the limit h — » 0. In the following 
subscripts o are used to indicate quantities which 
are are not yet our final estimators from the n th 
simulation. We define 

< +1 (i?) = e - S o +1 ^ C ^. 



H n (E) 



where the (otherwise irrelevant) constant c is in- 
troduced to ensure that Sq +1 (E) is an estimator 
of the microcanonical entropy 

S'£ +1 (E) = -\nc + S n (E)+\nH n (E) . (A.7) 

Inserting this relation into (A.3) gives 

(E) = (A.8) 



un+l 



b n (E) + [lnH n (E + e) - lnH n (E)]/e . 

The estimator of the variance of 6q +1 (_E) is ob- 
tained from 

a 2 [b n a +\E)]=a 2 [b n (E)} + 
a 2 [In H n (E + e)]/e + a 2 [In H n (E)]/e . 



Now a 2 [b n (E)} = as b n {E) is the fixed function 
used in the n th simulation and the fluctuations are 
governed by the sampled histogram H n — H n (E) 

<j 2 [\n{H n )] = [ln(H n + AH n ) - ln(iT)] 2 

where AH n is the fluctuation of the histogram, 
which is known to grow with the square root of the 
number of entries AH n ~ \J H n . Hence, under the 
assumption that autocorrelation times of neighbor- 
ing histogram entries are identical, the equation 



* 2 [b% +1 (E)] 



H n (E 



H n (E) 



(A.9) 



holds, where c' is an unknown constant. The as- 
sumption would be less strong if it were made for 
the energy-dependent acceptance rate histogram 
instead of the energy histogram. In the present 
models the energy dependence of the acceptance 
rate is rather smooth between nearest neighbors 
and there is less programming effort when using 
only energy histograms. Equation (A.9) shows that 
the variance is infinite when there is zero statistics 
for either histogram, H n (E) = or H n (E+e) = 0. 
The statistical weight for b^ +1 (E) is inversely pro- 
portional to its variance and the over-all constant 
is irrelevant. We define 



9 ° {E) ~ a 2 [b% +1 (E)] 



(A.10) 



H n (E + e) H n {E) 



H n (E 



H n (E) 



which is zero for H n (E+e) = or H n (E) = 0. The 
n th simulation is carried out using b n {E). It is now 
straightforward to combine 6q +1 (E) and b n (E) ac- 
cording to their respective statistical weights into 
the desired estimator: 

b n +\E) = g n {E)b n {E)+g%{E)b n Q +1 (E), (A.ll) 
where the normalized weights 



9o(E) 



and 



g*(E) + g$(E) 



g n (E) = l-g%(E) 



(A.12) 



(A.13) 



are determined by the recursion 

g n+1 (E)=g n (E)+gZ(E), g°(E) = 

We can eliminate b^ +1 (E) from equation (A.ll) by 
inserting its definition (A. 8) and get 



(A.14) 



b n+L (E) = b n {E)- 



(A.15) 



j%(E) x [\nH n (E + e) - In H n (E)]/e . 

Notice that it is now save to perform the limit ho — > 
0. Finally, equation (A. 15) can be converted into 
a direct recursion for ratios of the weight factor 
neighbors. We define 



R n {E) = e tbn{E) = 
and get 

R n+1 (E) = R n (E) 
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